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ABSTRACT 

Motivated by the observed connection between molecular hydrogen ( H2) and star formation, 
we present a method for tracking the non-equilibrium abundance and cooling processes of 
H2 and H2-based star formation in Smoothed Particle Hydrodynamic simulations. The local 
abundances of H2 are calculated by integrating over the hydrogen chemical network. This 
calculation includes the gas-phase and dust grain formation of H2, shielding of H2, and pho- 
todissociation of H2 by Lyman- Werner radiation from nearby stellar populations. Because 
this model does not assume equilibrium abundances, it is particularly well suited for simu- 
lations that model low-metallicity environments, such as dwarf galaxies and the early Uni- 
verse. We further introduce an explicit link between star formation and local H2 abundance. 
This link limits star formation to "star-forming regions," represented by areas with abundant 
H2. We use simulations of isolated disk galaxies to verify that the transition from atomic to 
molecular hydrogen occurs at realistic densities and surface densities. Using these same iso- 
lated galaxies, we establish that gas particles of 10* Mq or less are necessary to follow the 
molecular gas in this implementation. 

With this implementation, we determine the effect of H2 on star formation in a cosmo- 
logical simulation of a dwarf galaxy. This simulation is the first cosmological simulation with 
non-equilibrium H2 abundances to be integrated to a redshift of zero or to include efficient 
SN feedback. We analyze the amount and distribution of star formation in the galaxy using 
simulated observations of the HI gas and in various optical bands. From these simulated ob- 
servations, we find that our simulations are consistent with the observed Tully-Fisher, global 
Kennicutt-Schmidt, and resolved Kennicutt-Schmidt relations. We find that the inclusion of 
shielding of both the atomic and molecular hydrogen and, to a lesser extent, the additional 
cooling from H2 at temperatures between 200 and 5000 K increases the amount of cold gas 
in the galaxies. The changes to the ISM result in an increased amount of cold, dense gas in the 
disk of the galaxy and the formation of a clumpier interstellar media (ISM). The explicit link 
between star formation and H2 and the clumpier ISM results in a bluer galaxy with a greater 
spatial distribution of star formation at a redshift of zero . 

Key words: hydrodynamics, ISM: molecules, stars: formation, galaxies: dwarf, galaxies: 
evolution, galaxies: dwarf 



1 INTRODUCTION 

Strong observational evidence relates star formation to molecu- 
lar hydrogen, H2. For example, molecular clouds (MCs) are ob- 
served to be the sites of star formation both in the Milky Way (e. i 
Blitzlll993h an d in nearby galaxies jpukui & Kawamtiral I2OIC 
Wong & Bliti ('2002) showed that the star formation rate (SFR) 
in a set of CO-bright spirals is more closely correlated to the 
H2 abundance than the atomic hydrogen (HI) abundance while 



iFumagalh etal] ( |2009|) found that galaxies with low levels of star 
formation are observed to be H2 -poor. The connection between H2 
and star formation in other galaxies h as been further s upported by 



data from THINGS (Walter et al. 2008), HERACLES (Leroy et al 
^2009), NGS (Gil de Paz et al, .2007.). and SINGS (Kennicutt et al. 



I2OO3). Synthesizing these data. lLerov et all ( l2008h found that H2 
correlates strongly with the star formation efficiency (SEE). 

The best-known example of the connection between H2 and 
star formation comes from the Kennicutt-Schmidt (K-S) relation 
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( ISchmidl|[l959l : iKennicutll f 19981) . This law relates the surface den- 
sity of star formation to the surface density of HI, H2, or total gas 
and is frequently used as a benchmark for computational models 
of star formation. It has been calculated both over entir e galaxies 
and on smaller spatial scales. iBigieletalJ ( I2OO8II2OIOI) created a 
spatially resolved K-S relation using data from THINGS and re- 
lated surveys by finding the SFR surface density as a function of 
the gas surface density on sub-kpc scales. They found that at gas 
surface densities where H2 dominates (Egas > 9Mqpc~^) the 
SFE is higher than in lower surface density areas. This change in 
the efficiency implies that the SFR is directly related to the surface 
density of H2 and less directly to the total gas surface density. Ex- 
ten ding this study to gre ater radii in the galaxies by stacking spec- 
tra, jSc^^aetal] (|20Tl|) verified that this relationship held out to 
^25, where H2 had been previously undetected. Other studies of the 
K-S relation comparing the SFR to H2 at various spatial scales and 
across a wide sample o f galaxies also foun d that the SFR correlates 



most strongly with jRownd & Young| l999: Mur gia et al. 



2002 



2007 



Heyeretal. 2004; Gao & Solomon 2004; Kennic utt et al 
Blanc et al. 2009; Warren et al. 2010; Bigiel et al. 201^^ 

Theorists have developed different analytical models relating 
the SFR and H2 abundance to properties of the ISM. Based on 
the relationship between pressure and the transition fro m HI to 
H2 presented in lElmegreenl ( Il993h . IWong & BIit9 (|2002') hypoth- 
esized a dependency of H2 fra ction on midplane pre s sure for their 
sample of galaxies. Similarly, iBlitz & Rosolowskvl ( |2004 l2006h 
devel oped a pressure-based star formation law from this relation- 
ship. iKrumholz et al] j2012h argued for a volumetric, molecular 
K-S relation that would account for the variation in the total-gas 
surface density-based K- S relation observed betw een spiral, dwarf, 
and starbursting galaxies. IKrumholz et al.l ( l2009ah used an analytic 
model of the equilibriu m state of photodissociation fronts in MCs 
( IKrumholz et alj|2008l) to model the atomic to molecular transi- 
tion throughout galaxies as a function of metallicity and column 
density. This model was used to develop a star formation law ap- 
plicable to both primari l y atomic and primari l y mol ecular gas in 
IKrumholz et"aLl i2009bh . iMac Low & Gloved fcoiOl) . in contrast, 
presented a non-equilibrium, time-dependent analytic model for the 
H2 fraction in a region based primarily based on the metallicity and 
density. Analytic models of star-forming regions such as these have 
been extensively used to study H2 and star formation. However, 
simulations of galaxies are necessary to study how H2 evolves in 
concert with other processes. 

Until recently, simulations of galaxy formation have largely 
neglected H2 because of the high resolution required (see below 
for a description of more recent simulations that do include H2). 
Including H2 in galaxy formation simulations, however, has the 
potential to affect both the SFR and the spatial distribution of star 
formation. If the formation of H2 is a fundamental step in star for- 
mation, then this omission would result in a less physical connec- 
tion between the gas properties and star formation. Linking star for- 
mation to H2 would increase the connection between star forma- 
tion and the properties of the gas that determine the fraction of H2. 
In addition to volume density, these factors may include the gas 
surface density, the metallicity, and the dissociating radiation field. 
The distribution of both metallicity and the radiation field could af- 
fect the distribution of star formation and should potentially reflect 
this variation. As will be discussed further, reproducing the coiTect 
spatial distribution of stars is important, not only for comparisons to 
observations, but because stars can, through supernova (SN) feed- 
back, impact both the baryonic and dark m atter (DM) mass dis- 
tribution of galaxies jGovemato et al.l [20101) . Furthermore, these 



properties vary across galactic type and redshift, potentially lead- 
ing to varying SFRs. For this reason, an H2-based star formation 
law has been used to motivate reduced star formation efficiencies in 
the low-metallicity en vironments of dw a rf galaxies and at high red- 
shift (e.g. [R obertson &Kravtsovl [20081 ; iGnedin & Kravtsovll20ld ; 
lKrumholz & Dekel 2oTTx 

A second potential effect of including H2 is in the modeling of 
the cold, dense ISM. The gas shielding necessary for H2 formation 
reduces the photoheating, leading to increased amounts of cold gas. 
The importance of shi elding to the creation o f col d, star forming- 

i ;as ha s been stressed in lKrumholz et al.l l l201lh and lGlover & ClarM 
201 ih . Furthermore, while H2 does not em it strongly at the cool 
temperatures of MCs jGlover & Clarkl201 l|), H2 can in low metal- 
licity regions be an important coolant at temperatures between 
5000 and 200 K through coUisions ( Glover & Mac Low 20(33; 
[Onedin & Kravtsovlboioh . In the early universe this c ooling was 
respo nsible for the formation of Population III stars dAbel et al.l 
I1997I) . At early times and in low-metallicity environments, it could 
also have a strong impact on the density structure of the gas and 
the am ount of gas and star formation in DM halos ( Jappse n"et al.l 
l2007h . 

The combination of shielded gas and higher low-temperature 
cooling rates promotes the formation of cold, dense clouds. Includ- 
ing H2 is, therefore, important for better reproducing the environ- 
ment where star formation is observed to take place. The forma- 
tion of a clumpy ISM has been shown to increase the efficiency of 
SN feedback at removing gas, jGovemato et al.ll2010l ; lBrook et al.l 
I2OI0I I2OIII) . Consequently, by enabling the formation of cold, 
dense gas and limiting star formation to these star formation re- 
gions, H2 has the potential to make the ISM more prone to out- 
flows. Outflows such as these can change both the mass and mass 
distribution of galaxies. Because of this interaction between feed- 
back and the thermodynamics of the ISM, it is necessary to include 
both H2 and efficient SN feedback in simulations. 

In the past five years, several models have been developed 
to eith e r stud y MCs or to model H2 i n simu l ations of galaxies. 
iDobbsl j2008h ; iDobbs et alj ( I2OIII) and lTaskeJ ( l201lb studied the 
formation, orbits, and mass distributions of MCs in models of spi- 
ral galaxies by identifying gas clumps where MCs would be. These 
simulations followed the dynamics of MCs but do not include the 
formation or destruction of hydrogen molecules and the heating and 
cooling associated with them. 

Other simulations of galaxy evolution have included models 
for the H2 abundance and used i t to calculate the SFR. Based 
on observed scaling laws for MCs ( lLarson|[l98lb . IPelupessv et al.l 
( 2006) presented a time-dependent sub-grid method for simulat- 
ing MCs in Smoothed Particle Hydrodynamic (SPH) simulations. 
This method has been used to reproduce the effect of metal- 
licity on the H2 abundance in dwarf galaxies dPelupessv et al.l 
l2006h . the relationship between H2 abundance and mid-plane 
as pressure, and the con versio n factor between CO and H2 
Pelupessv & Papadopouloj|2009l ). Because H2 was linked to star 
formation, this model has also been used to study deviations 
from the K-S relation in low-me t allicity, high-redshift galaxies 
dPapadopoulos & Pelupess\ll201(]|) . iRobertson & Kravtsovl (2008) 
also tracked the local H2 abundance in their simulations of galax- 
ies. For their model, they used the equilibrium abundance of H2 
for gas particles, based on the metallicity of the gas and the amount 
of dissociating radiation. With these abundances, they were able 
to examine the importance of H2 to star foiTnation in their SPH 
simulations of isolated galaxies and reproduce deviations from the 
global K-S relation in low metallicity galaxies based on reduced 
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H2 fractions. i Kuhlen et a l] llOll) applied the equilibrium mod- 
els for H2 abundance from Krumholz et al. ( 2009a) to models of 
dwarf galaxies at high z and found that in the absence of efficient 
SN feedback, H2-based star formation reduced the stellar mass in 
dwarf galaxies. Each of these models, however, rely either on as- 
suming the global properties of MCs or assuming that the gas is in 
chemical equilibrium. 

In the o nly f ully self-consistent implementation so far, 
iGnedin et alj ( l2009h implemented a model for the local non- 
equilibrium H2 abundances in their Adaptive Mesh Refinement 
(AMR) simulations. In this model, they used local rates of H2 for- 
mation and destruction to integrate over the chemical network. Us- 
ing this method, they created cosmological simulations of galaxies 
to determine deviations from the K-S relation in low metallicity 
galaxies (Gnedin & Kravtsov 2010) and with different dust-to-gas 
ratios and far-UV fluxes (Gnedin & Kravtsov! l201lh and to study 
sources of scatter in the K-S relation jFeldmann et al fcoilh . This 
model has since been used to verify equilibrium analytic mod- 
els for the H2 abundance for gas with metallicity > 10"'^ Zq 
jKrumholz AGn edin 2011). 

In this paper, we have adapted the method presented in 
iGnedin et alj l l2009l) for SPH simulations with efficient SN feed- 
back. We present a method for integrating the local non-equilibrium 
abundance of H2 in SPH simulations based on the local formation 
and destruction rates. These rates include the formation of H2 on 
dust grains and photodissociation of H2 by Lyman Werner (LW) 
radiation from nearby stellar populations. The calculation of the lo- 
cal LW flux is based on the radiation from nearby stellar particles, 
in which the gravity tree structure is used to establish proximity. 
Accordingly, we also model the shielding of H2 from the LW radia- 
tion by both itself and dust for each gas particle. We also include the 
dust shielding of atomic hydrogen from the cosmological UV back- 
ground. Molecular hydrogen-based star formation and a blastwave 
model for SN feedback are also included in the simulations, in or- 
der to produce both star formation and SN feedback in regions with 
high H2 abundances. This model is ideally suited for simulations 
of galaxy formation because of the speed and ability to simulate 
large spatial ranges inherent in SPH. By integrating over the chem- 
ical network, no assumptions are made about the final structure of 
molecular clouds or the relationship between H2 and global galac- 
tic properties. It is, therefore, especially useful for studying extreme 
environments, such as the early Universe and low-metallicity dwarf 
galaxies. 

Despite the work done on H2 in evolving galaxies, few of the 
previous simulations with H2 have been cosmological and none 
have been integrated to a redshift of zero. Furthermore, none of 
them have included both H2 and efficient feedback. Cosmologi- 
cal simulations are necessary to determine the effect of H2 on star 
formation and disk structure in realistic models of galaxy forma- 
tion that include gas accretion through both cold-flows and merg- 
ers. Following these simulations to a redshift of zero is necessary 
to analyze the effect of H2 on the stellar and gaseous content of 
galaxies over the history of the Universe and to compare them to 
observations of nearby galaxies. The inclusion of efficient feedback 
is important for correctly modeling star formation and gas loss and 
for examining the connection between the formation of cold, dense 
gas and outflows. For these reasons, we present a cosmological sim- 
ulation of a dwarf galaxy with non-equilibrium abundances of H2 
and efficient feedback integrated to z=0. 

We focus this work on field dwarf galaxies, both because of 
their manageable computational scale and scientific interest. Dwarf 
galaxies are well suited for simulations because they require fewer 



particles to achieve high mass resolution. Their low virial masses 
and metallicities also make them sensitive laboratories for study- 
ing gas dynamics and star formation. Having a small virial mass 
means that gas heated through SN feedback can more easily escape 
the halo. Their low metallicities reduce the cooling rates of the gas, 
making H2 a relatively more important coolant. Their relative lack 
of dense gas also makes these galaxies an extreme environment for 
star formation. Their SFE may be even further reduced by the dif- 
ficulty of H2 formation in such a low metallicity and low surface 
density environment. With their low metallicities, the formation of 
H2 and, consequentially, star formation will occur at higher densi- 
ties than in spiral galaxies. 

This method may be modified for other SPH simulations and 
we anticipate that our results will be of interest to those concerned 
with modeling the connection between ISM and star formation and 
to those interested in star formation in dwarf galaxies. In !j2l we de- 
scribe our algorithms for calculating H2 abundance and simulating 
star formation and test our model on isolated disk galaxies. In §3 
we present the results comparing the effect of different models of 
the ISM and star formation recipes on the amount and distribution 
of star formation. Finally, in §4, we discuss how the changes in the 
star formation were related to changes in the thermodynamics of 
the gas and discuss the connection to other works. 



2 METHODS: MOLECULAR HYDROGEN 
IMPLEMENTATION 

We implement a method for tracking the local, non-equilibrium 
abundance of H2 in our SPH simulations based on the local for- 
mation and destruction rates. The abundance of H2 is dominated 
by the balance between dust grain formation (§ 12.2.1b and pho- 
todissociation rates. We estimate the photodissociation rate based 
on the LW flux from surrounding star particles (i; 12. 2. 2b . The gas 
is in turn shielded from LW radiation by itself and dust (ij I2.2.3I >. 
Gas-phase reactions (i) l2.2.4b also impact the abundance of H2 in 
low-metallicity environments and provide an additional source of 
cooling. We model each of these processes and use a semi-implicit 
solver to calculate the H2 abundance in each gas particle at each 
time-step. This impleme ntation for ca l culat ing H2 abundances is 
based on the approach of Gnedi n et al.l ( 12009 ) but with a simplified 
estimate for LW radiation and applied to SPH, rather than AMR, 
simulations. These simulations also include a more efficient SN 
feedback model. It is the first self-consistent non-equilibrium H2 
implementation to be used in conjunction with efficient SN feed- 
back and in cosmological simulations integrated to 2 = 0. 

We tested our implementation on isolated models of a Milky 
Way-mass galaxy (i; 12.31 and ij |2.4b before completing the cosmo- 
logical dwarf galaxies. We, therefore, assess our H2 implementa- 
tion over a range of masses in this paper. The initial conditions for 
the isolated simulations are described in Appendix A. In particu- 
lar, we ensured that the transition from atomic to molecular hydro- 
gen occurred at the correct densities and surface densities at our 
resolution and displayed a realistic metallicity dependency. This 
transition is an important constraint because the presence of H2 is 
observed to be strongly correlated with the surface density of gas. 
Furthermore, many theoretical models for the abundance of H2 pri- 
marily depend upon the extent that the gas is shielded, a quantity 
that is directly related to surface density. We used these same sim- 
ulations to determine the necessaiy resolution to model H2 and the 
conditions in which it forms. 
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2.1 General Description of Code 



where Xi is the number fraction of baryons in species i 



We implement H2 physics in the galaxy formation code , GASO- 
LINE. GASOLINE is a parallel tree+SPH code ( Wads lev et all 
[2004) and an extension of the A'-body gravity code PKDGRAV 
( IStadelll2001l) . GASOLINE has been used extensively for studying 
galaxy formation and evolution an d has been successful in repro- 
ducing the TuUy-Fisher Relation (Govemato et al. 2009), the ro- 
tation curves of dwarf galaxie s (Go vemato et al. 2010) , the inter- 
nal st r ucture of disk galax i es jGovem ato et al. I2007t iRoskar et al.l 
l2008l : IStinson et al] l2009l : iBrooks et al., ,20091) . the column den- 
sity distrib utions of damped Lyman- a systems jPontzen et al.l 
I2OO8', 2010), and the observed stellar mass-metallicity relationship 
iBrooks et al. 2007). 

In addition to gravity and gas dynamics, GASOLINE models 
physical processes such as ion abundances, background UV, gas 
cooling rates, star formation, and feedback from SNe. We use a 
semi-implicit integrator to calculate the non-equilibriu m ion abun- 
dance s and gas cooling from coUis ional ionizati o n rates ( Abel et al 
19971) , radiative recombination jBlackl Il98lt IVemer & Ferlanc 
\99&), photoionization, bremsst rahlung, and cool ing from H, He 
(Cen 1992) and metal lines jShen et al.l I2OI0I) . For the pur- 
poses of calculating the rates of radiative recombination and pho- 
toionization, our simulations include a redshift-dependent UV- 
background radiation fr om quasars and stars in external galaxies 
jHaardt & Madatjl 19961) . Cooling rates from metal lines are tabu- 
lated based on the gas temperature, density, and metallicity and the 
UV backgr ound. These rates we re computed using CLOUDY (ver- 
sion 07.02; iFerland et al.|[T998h with the assumptions that the gas 
is in ionization equilibrium and optically thin to UV radiation. In 
this work, we extended these calculations to include the presence 
of H2, as described in i]2.2| 

Star formation in GASOLINE is determined on a probabilistic 
basis accordin g to the free-fal l time of potentially star-spawning gas 
particles (Stin son et al.ll2003) . We introduce a modification to our 
star formation recipe, further described in ^2.5\ in which the abun- 
dance of H2 is incorporated into the star formation probability. In 
order to calculate the effects of feedback fro m Type la and Type I I 
SNe, we use the blastwave method outlined in I Stinsonetal1j2006l) . 
In this method, energy from SNe is distributed among surrounding 
gas particles and gas cooling is disabled for a length of time and dis- 
tance determined by the blastwave model. This feedback method is 
very efficient at producing gas outflows. Supemovae are also re- 
sponsible for enriching the surrounding gas particles and metals 
are distributed th roughout the simulation through metal diffusion 
IShen etal.ll2010l) . 



2.2 Implementation of Molecular Hydrogen 

The abundance of H2 is based on the rates of dust grain H2 for- 
mation, photodissociation with shielding, gas-phase H2 formation, 
and collisional dissociation. The abundance of HI is now a function 
of the rate of H2 formation, in addition to the rates of photoioniza- 
tion, recombination, and collisional dissociation. Using these rates, 
the equations of balance for HI and H2 may be written as 



Xui = R{T)neXuii — SdXninijTui — Chit^hiWe — 2^h2 

(1) 

~ RdnbX-H_i{X-Hi + 2XH2) " SdSn^X u^nbY^^ + ^ifj 

(2) 



Xi — Hi /rib, 



(3) 



R{T) is the recombination rate, rie is the number density of elec- 
trons, Sd is shielding from dust, nj, is the number density of 
baiyons, Fhi is the hydrogen photo-ionization rate, Chi is the HI 
collisional ionization rate, Rd is the formation rate coefficient of 
H2 on dust, S'h2 is self-shielding of H2, is the H2 photodis- 
sociation rate, and X^^^ is the gas-phase dissociation and formation 
rate of H2. The shielding factors, Sd and Sn^, represent the frac- 
tion of ionizing or dissociation flux that reaches the gas and range 
from zero to one. The gas-phase dissociation and formation rate in- 
cludes collisional dissociation of H2 and the gas-phase formation 
of H2 via H~. We describe the dust-based formation, the rate of 
photodissociation and the shielding, and the gas-phase dissociation 
and formation in the following sections. 

2.2.1 Formation on Dust Grains 

The formation of H2 on dust grains dominates over gas-phase 
formation once a s mall reservoir of metals is available. As in 
iGnedin et alj ( l2009h . we use an observationally motivated equa- 
tion for the coefficie nt for rate of H2 formation on dust, Rd, from 
IWolfire et al.l ( I2OO8I) . In this formalization, the rate of H2 forma- 
tion on dust depends on both the amount of dust and the gas den- 
sity structure. In order to calculate the amount of dust, we assume 
a dust-to-gas ratio proportional to the gas metallicity, Z. 

The sub-grid gas density structure affects the rate of H2 
formation through the increased rate of formation in unresolved 
high-density regions. Therefore, any model that follows the non- 
equilibrium abundances of H2 without resolving the turbulent 
structure within gas clouds must somehow parametrize the struc- 
ture. W e address the sub-grid nature of H2 formation in the same 
way as iGnedin et al.l (2009), by assuming that the gas has a log- 
normal density substructure. On scales smaller than gas particles, 
the increased formation rate caused by the substructure can be char- 
acterized by a clumping factor, Cp, equal to (p^) / (p). For a log- 

2 

normal density distribution Cp = e '"p, in which Up is the width 
of the distribution. Within our simulation s Cp is an adj ustable pa- 
rameter with a best-fit value of lO.O, as in Gnedin et al. (2009). 
The accuracy of estimating the formation enha ncement with 



this c lumping factor has been confirmed by iMilosavlievic et al 



2017) for H2 fractions up to 50%. However, Milosavlievic et al 



2011) also cautioned it may over-predict the rate of formation at 
higher densities and H2 fractions. Therefore, further development 
to the equation for H2 formation at high densities and molecular 
fractions may be waiTanted. However, we choose to focus our anal- 
ysis on situations in which fully molecular gas is rare (dwarf galax- 
ies) and where the transition from mostly atomic to predominantly 
molecular gas is more important than the exact H2 abundance at 
very high densities. Including both the effect of dust and sub-grid 
clumping, our H2 formation rate is defined as 

3.5 X W'^'^Z/ZqCp cm^s"^ 



Rd 



(4) 



2.2.2 Photodissociation by Lyman Werner Radiation 



Dissociation of H2 is caused primarily by LW radiation (F jj^ )^ 
which can either be produced by young stars in the galaxy or can 
come, to a much smaller degree, from extra-galactic sources. We 
assume the LW radiation from extra-galactic sources takes the form 
of a cosmological background varying with redshift and use the 
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function from lHaardt & Madaul ( Il996 l). At 2 = 0, this radiation is 
approximately 5000 times less than the average flux of LW radia- 
tion in the disks of galaxies. For z > 2, the background radiation 
approaches 50 times less than the average flux. We include cosmic 
LW radiation for consistency but it has a negligi ble impact on the 
star formation, as also observed in GnedinldlOlol) . 

Stellar light in the galaxy is the primary source of LW radia- 
tion. This radiation depends on the distribution of young stars and 
can only be correctly determined using full radiative transfer. Un- 
fortunately, full radiative transfer is currently computationally pro- 
hibitive in cosmological simulations that are evolved for a Hubble 
time. Even those cosmological simulations that do include radia- 
tive transfer must make assumptions, such as using an optically thin 
approximation to find the Eddington tensor dOnedin & Abelll200ll : 
[Petkova & Springel 2009). 

For the sake of speed and computational ease, we make a sim- 
ple approximation of the amount of stellar LW flux a given gas 
particle experiences based on the average flux from nearby stars. 
When calculating the flux, we assume a completely optically thin 
disk. However, as will be shown in i) l2.2.3l we do include shielding 
from both dust and H2 molecules when calculating the photodis- 
sociation rate of H2 for a gas particle. This method for calculating 
the LW flux relies on LW dissociation being primarily important 
over short distances and the absorption of LW photons occurring 
mostly in MCs. Fortunately, the efficient shielding of H2 at high 
densities causes the fraction of H2 to be relatively insensitive to 
the amount of dissociation rad i ation incident on MC s dOnedin et al.l 
I2OO9I : iKmmholz et alj|20093 : iMac Low & Gloved l2010h . As de- 
scribed below, this method is sufficient to follow the large-scale 
spatial and temporal variations of the LW radiation and to repro- 
duce the surface density at which hydrogen becomes molecular to 
within a factor of two. 

The details of our method for calculating the LW flux are as 
follows. The average flux impacting a gas particle from nearby stars 
is calculated using the tree built for the gravity calculation. To do 
this, we apply the flux from star particles within the same cell or, 
if necessary, the same parent cell in the following manner. Within 
each cell of the tree, we assume that the amount of LW flux coin- 
cident on a given gas particle is a function of the amount of LW 
radiation emitted by star particles within the same cell. Each star 
particle emits LW radiation ( L) according to its a^ e (f) and mass 
(m), based on the Starburst99 (Leitherer et al. 1999) calculation for 
a simple stellar population with a Kroupa et al. (1993) IMF. We 
then define one average LW source for each cell with luminosity, 
Ls, equal to the total luminosity of stars within the cell and located 
at the luminosity-weighted average location, Xs, of the star parti- 
cles. Based on the luminosities of individual star particles and their 
locations (xi), we define Ls and Xs as 

N 

Ls^^L{U,mi) (5) 

Xs — \p) 

Lis 

Gas particles within the cell experience a LW flux as if all ra- 
diation were coming from the average LW source. For a gas particle 
with a distance greater than one softening-length away from the av- 
erage source, we set the flux equal to what it would be from a point 
source at that distance and with that luminosity. For a gas particle 
with a distance less than the gravitational softening, we set the flux 



equal to the flux at the softening radii. The flux, F, incident on a 
gas particle at position Xg is, therefore, 

I . I — n iflxs— Xgl>/i 
1 . f -I otherwise 

where h is the softening. This adjustment prevents artificially high 
levels of radiation at distances too small to be resolved. 

One complication to this method arises when there is a bright 
LW source nearby but not within the cell. If the radiation from this 
source were ignored, the amount of LW radiation on the gas par- 
ticles in the cell would be underestimated. In order to account for 
cases in which bright sources in adjacent cells should impact the 
amount of flux, we make the following adjustment. We use the av- 
erage LW source of the parent cell if doing so results in the gas 
generally experiencing higher amounts of LW radiation. 

In order to estimate whether the radiation would be generally 
higher using the parent cell average source, we calculate the loca- 
tion and luminosity of the average sources for both the parent and 
child cell. We then determine what the flux from the parent average 
source would be at the center of mass of the gas in the child cell. 
We compare this to the flux from the child average source at a dis- 
tance equal to the first moment of the gas mass in the child cell. If 
the flux is higher in the first case, we calculate the amount of radia- 
tion experienced by each gas particle in the child cell based on the 
average LW source of the parent cell. Otherwise, we only use the 
data for the child cell and proceed as in equation 7. 

We illustrate the accuracy of our method using a simulation 
of an isolated Milky Way-mass galaxy (MWmr from Appendix 
A). Figure [T] compares the LW flux at each gas particle calculated 
within the code to the theoretical LW flux, assuming an optically 
thin ISM.Q In our simulations, this method produced fluxes within 
a factor of ten of the optically thin result. This factor of ten equates 
to a change of about 50% to the density and surface density at 
which the transition from HI to H2 takes place. Despite the lim- 
itations of our method, it is able to track general trends in the LW 
radiation caused by different rates of star formation. 

2.2.3 Shielding from Radiation 

Shielding from LW radiation is the primary reason molecu- 
lar gas can form and remain in high-density regions. We used 
a phenomenological model to approximate the computation- 
ally expensive three-dimensional treatment of shielding based 
on the surface density of each gas particle, following the 
work oflDra ine & Be rtoldH ( ll996l) , lGlover & Mac Lowl ( |2007|) and 
iGnedin eta l. (2009). To fully reproduce the shielding would re- 
quire a three-dimensional treatment of radiative transfer and is be- 
yond the computational capabilities of our cosmological simula- 
tions. However, this simplified model is able to reproduce the ob- 
served efficiency of shielding from LW radiation as a function of 
surface density and metallicity. 

In this model, the amount the gas is shielded from LW radia- 
tion is a function of the gas surface density and metallicity (used to 

^ To calculate the theoretical amount of LW flux for comparison, we first 
calculated the amount of LW each stellai' particle emitted based on the Star- 
burst99 (Leitherer et al. 1999) value for simple stellar population with the 
same initial mass function (IMF), mass, and age. Then, for each gas par- 
ticle, we summed the amount of flux from each stellar particle at the gas 
particle's location, assuming there was no intervening absorption. 
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Figure 1. Comparison between our approximation of the LW flux experi- 
enced by eacli gas particle in simulation MWmr (Code FIuxlw) and the 
amount of flux each gas particle would experience in the optically thin limit 
(Calculated FIuxlw)- The histogram is normalized by the total number of 
particles. 



calculate the dust shielding) and H2 surface density (used to cal- 
culate the self-shielding). We also apply dust shielding to HI when 
calculating the rates of photoionization and heating from the UV 
background radiation. The functional forms for the dust shielding, 
Sd, and self-shielding, S'hj. are based on observations and the pa- 
rameters were tuned to fit observational data. 



= 



1 — CJH2 , -0.00085(1+1)1/2 



+ 



+ (l + a;)(i/2) 



(8) 
(9) 



In the above equations, A^hi is the column density of HI, is 
the column density of H2, and x = N-h2/{5 x 10^'*) cm^. Both 
Cd^ff and UJ-H2 are adjustable parameters tuned to be 4 x 10^^^ cm^ 
and 0.2, respectively. 

The HI and H2 column densities are calculated based on the 
volumetric number density of a gas particle, rib, and an assumed 
column length of the particle. We set the column length for both 
dust and self-shielding to be equal to the smoothing length, h, of 
the gas particle, resulting in Ni = Xi h rib, where i may be either 
HI or H2. Our use of h for the column length was largely decided 
on for empirical reasons (see § 12.31 for a test of our H2 calcula- 
tion, including this assumption of column length). However, it also 
has the advantages of simplicity and being directly related to the 
volume over which the gas particle mass is spread. We also inves- 
tigated the use of more physically motivated equations of the col- 
umn length, including the use of a turbulent sc ale length based on 
the local shear, as presented in IPavlovski et al.i (2Q02). We found, 
however, that in SPH simulations, basing the column length of a 
gas particle on the shear across a particle resulted in unphysically 
small column lengths at the densities where H2 becomes abundant. 



2.2.4 Gas-Phase Dissociation and Formation 

Gas-phase reactions, including both coUisional formation and dis- 
sociation, can be important for establishing the abundance of H2 in 
low-metallicity environments and can provide an additional source 
of cooling for the gas. Gas-phase formation of H2 occurs via 
+ H ^ H2 + H+ an d H~ + H -> H2 + e". Following the 
work of I Abel et al. I ( Il997t) , we assume that the abundance of H2 



is largely independent from the abundance of H and use equation 
26 from their minimal model for H2 formation via H^. In addi- 
tion to photodissociation, H2 may be dissociated through collisions 
with H2, HI, HII, and e^. We use the rates for these collisions 



given in ITepp & Shull ( 19831) , lOove & Mandvl ( ll98d) . lAbel et alj 
(l99f), and Do nahue & Shulf jl99lh . respectively. These processes 
together result in 

^ kin^-nu—nn2ik2nu2+k3nH + k4nnu + k5ne) (10) 

where rijj- is as defined in equation 27 of I Abel et al. I ( Il997h . ki 
is the H2 formation rate via njj- , and ^2, fca, ^4, and fcs are the 
coUisional dissociation rates by H2, H, HII, and electro ns, respec- 
tively. The values for ka, k^, and fcs are summarized in lAbel et al.l 
( I1997I) as reactions 11, 12, and 13, respectively. For all gas-phase 
calculations, we assume a constant ratio of ortho-hydrogen (nuclear 
spin number / = 0) to para-hydrogen (J = 1 ) equal to 3:1. This 
assumption is justified in lGlover &Abell ( l2008h . 



2.2.5 Heating and Cooling 

With the implementation of H2, we added cooling of H2 due 
to coUisional dissociation and collisionally excited line emis- 
sion. CoUisional reactions are particularly important for the cool- 
ing of H2 when 200 K < T < 5000 K and p > 10 amu/cc) 
jGnedin & Kravtsovll2010l) . 

Dissociative cooling due to collisions with H2, H, HII, and 
electrons is calculated usin g the rates and amount of emission pre- 
sented in lAbel et al. I( ll997l) . The coUisional induced line-excitation 
cooling rates are determined for interactions between H2 and H2, 
H, HII, Hel, and electr ons. For these interactions , we use the 
cooling rates presented in iGlover & Mac Lowl ( I2OO7I) for the low- 
density case and summarized in Tabl e 8. The line-excitation cool- 
ing rates from Glover & Mac Low ( 2007h use value s calculated in 
IWrathmall et al.l ( l2007i) from the iMielke et al.l ( l2002l) potential sur- 
face an d are substantially gr eater than the previously determined 
rates of iGaUife Fallal ( 1 19981) . Deuterium cooling is ignored as it 
is primarily important at densities too high to be resolved in our 
simulations. In addition to photoionization heating for atomic and 
ionized hydrogen and helium, we include heating due to photodis- 
sociation of H2. 



2.3 Testing Molecular Hydrogen Abundances 

In order to verify the accuracy of our H2 abundance methodology, 
in particular the shielding calculation, we compared the volume 
densities and surface densities at which gas transitioned between 
atomic and molecular in our isolated simulations of different metal- 
licities. We found that the transition occurred at densities between 
p = 10 amu/cc and 100 amu/cc, depending on the metallicity, 
as predicted (Figure |2). These values are consistent with average 
densities of MCs. 

The surface density at which the transiti on occurs is shown in 
Fig ure [3] and compared to observations from lGillmon et all (|2006|) 
and lWolfire et al.N2008h . When determining the dependence of the 
H2 abundance on surface density, we show the surface density used 
in the code calculation (Ni — Xi h nt) (grey contours) and, in 
the case of MWmr, the surface density as it would be determined 
observational l y by the method used in iGillmon et al. I JioOQ) and 
IWolfire et al.l J2OO8I) (black contours). In these studies, the H2 sur- 
face density was determined from the amount of absorption in the 
quasar spectra. The HI surface density was determined by match- 
ing 21 cm emission with the same velocity line profile as the H2 
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Zg 0.5 Ze 0.1 Ze 




Density [omu/cc] 

Figure 2. The transition from atomic to molecular hydrogen for three iso- 
lated, Milky Way-like simulations with different metallicities. The points 
represent individual SPH gas particles. Vertical lines at 100 amu/cc (an av- 
erage MC density and the minimum density for star formation in our non- 
H2 simulations) are included to guide the eye. As the metallicity decreases, 
the combination of longer formation times and less shielding causes H2 to 
form at increasingly high densities. 

gas. In order to mimic this method, we post-processed the galaxy 
to create velocity cubes of the H2 and HI gas with a 2.5 pc by 2.5 
pc by 10 km/s grid. The velocity resolution was chosen to contain 
the typical line width of a MC. For each cube in the grid, we then 
determined X112 and A'^hi + '2Nu2- The correspondence between 
our determination of the surface density and the observational data, 
particularly for large fractions of H2, supports our use of h for the 
column length. 



Zo 0.3 Ze 0.1 Ze 




l°g(SHi + H.) [amu/cm^] 



Figure 3. Transition from atomic to molecular hydrogen as a function of 
surface density for three isolated. Milky Way-like simulations with dif- 
ferent metallicities. All th e contours rep r esent data from the simulations. 
The crosses are data from iGillmon et al.l | |2006|) and the X's are data from 
IWolfire et al ] j2008l) . The grey contours were calculated assuming surface 
densities of Ni = Xi h p, as used in the code. They represent the num- 
ber of points (individual gas particles) enclosed and are linearly spaced. 
The black contours in the left plot were calculated to mimic observations of 
MCs. For these contours, the surface densities were calculated for each cell 
in a 2.5 pc X 2.5 pc X 10 km/s velocity cube, in which 10 km/s was cho- 
sen to be an average velocity dispersion of a MC. These contours represent 
the distribution of the individual cells in surface density and H2 -fraction 
space and are linearly spaced at the same levels as the grey contours. The 
transition between atomic and molecular hydrogen in the solar metallic- 
ity simulation fits with the observational data. As metallicity decreases, the 
transition moves to higher column densities as both dust shielding and the 
rate of H2 formation on dust decrease. 



2.4 Resolution 

Over what resolutions does this implementation of H2 hold and 
over what resolutions are high-density regions sufficiently well re- 
solved to use it? When simulating H2, the mass and spatial resolu- 
tion has the potential to affect both the densities of the gas particles 
and their smoothing lengths, h. Molecular clouds have typical den- 
sities of 100 amu/cc, so the density structure within the simulations 
must be sufficiently resolved for gas with such densities to exist. 
The resolution of the density structure is a function of both the 
mass resolution (number of particles) and the force resolution (the 
softening length and smoothing length). 

Since the column densities used when calculating shielding 
are a function of the SPH smoothing, h, any resolution effect on 
them can also affect the fraction of H2. The value of h is primarily 
a function of the particle number but it has a minimum value of 
0. 1 times the softening length. In all of our simulations, we ensured 
that the minimum smoothing length was small enough that giant 
MCs could be resolved. 

In order to determine the effect of mass resolution on the 
amount of H2 formed, we simulated the same disk galaxy with 
lO'' (MWlr), 10*^ (MWmr) and 10^ (MWhr) gas particles. The ini- 
tial gas particle masses for each of the simulations were, therefore, 
of 1.4 X 10*^, 1.4 X 10*^ and 1.4 x 10* Mq. These mass resolutions 
were previously found to be sufficient to model the star formation 
and density structure in spiral galaxies prior to the addition of H2 
dChristensen et al.l20IClh . The effect of mass resolution on the frac- 
tion of H2 as a function of density is shown in Figure |4] We found 
that the relationship between density and the H2 abundance was 



independent of the number of particles in the three simulations, in- 
dicating that the transition from HI to H2 is consistent over this 
range of mass resolution. This figure also shows that for our range 
of resolutions, the density structure is sufficiently resolved for gas 
with densities up to 1000 amu/cc to form. 

However, the resolution does affect the fraction of gas that 
reaches high density (Figure |5)- As the resolution increases, the 
density distribution of gas particles becomes a smooth, declining 
function of density. We, therefore, recommend particle masses of 
1.4 X IC* Mq or less when using this method to simulate H2. 



2.5 Star Formation 

In our simulations with H2 we incorporated a link between star 
formation and MCs by modifying the probability of a star form- 
ing to include the local fraction of H2. This link limits star forma- 
tion to "star-forming regions," represented by areas with abundant 
H2. We base d our star formation recipe on the formalization of 
IStinson et al] ( 12006) . In this recipe star formation occurs stochas- 
tically, based on the local K-S relation in which the probability of 
star formation is a function of the local density. Gas particles that 
are sufficiently dense (p > pmin) and cool (T < Tmax) have a 
probability, p, of spawning a star particle. The value of p is a func- 
tion of the local dynamical time tform'- 

p= ^^(1-6-'=*^*''*/— ) (11) 
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Figure 4. The transition from atomic to molecular hydrogen for three iso- 
lated, Milky Way-like simulations with different mass resolutions. The 
points represent individual SPH gas particles. For each of the simulations, 
the transition occurs at the same density, indicating that within this range of 
resolution the transition is independent of the number of particles. 




log p [amu/cc] 



Figure 5. Normalized histograms of the gas densities within the three iso- 
lated. Milky Way-like simulations with different mass resolutions. 

where rUgas is the mass of the gas particle, mstar is the initial mass 
of the potential star particle, and c* is a star forming efficiency 
factor. 

In order to incorporate H2, we modified c* such that 

c* = Co (12) 

With this modification, the probability of star formation increases 
in highly molecular regions. In our simulation without H2,c* — Cq 
= 0.1, whereas with H2, Cq =0.1. 

In addition to incorporating into the star formation effi- 
ciency, we also altered the qualifications for a gas particle to form 
stars. In our non- H2 star formation r CCipC, pmin = 100 amu/cc, 
in order to mimic star formation in MCs. In our H2 star formation 
recipe, the dependence of c* on the molecular fraction already lim- 
its star formation to dense, molecular gas. We, therefore, reduced 
pmin to 0. 1 amu/cc, which is sufficiently low to have negligible ef- 
fect on the star formation. Similarly, we also lowered Tmax ■ In sim- 
ulations without H2, low-metallicity gas cannot cool below 10000 
K so we set Tmax ~ 20000 in order to select for the gas that would 
be in the warm and cold ISM. When H2 is included, however, a 



cold phase of the ISM is able to form. In order to select for gas that 
had cooled, we set Tmax = 1000 but found the star formation to be 
largely insensitive to this value. The values for these star formation 
parameters are also outlined in Table 2. 

When calculating star formation, the ability o f the simulation 
to reso lve the Jeans length and mass can be a factor. Isate & Burkerll 
( I1997I) showed that when the Jeans mass or and length of a parti- 
cle exceeds the particle mass and smoothing or softening length, 
respectively, star formation becomes dependent on the specifics of 
the star formation recipe. We chose not to include an artificial pres- 
sure floor to support th e gas against gravitational col lapse as other 
authors have done (e.g. iRobertson & KravtsovlboOSi) . Our simula- 
tions are sufficiently high resolution that gas particles for which the 
Jeans mass or length are unresolved necessarily meet our star for- 
mation criteria (including H2 abundance). Given Cq = 0.1, these 
particles will soon form stars, regardless of whether or not a pres- 
sure floor is included. Therefore, Jeans instability below the resolu- 
tion of our simulations is included within our sub-grid star forma- 
tion model. We found in several test simulations that the inclusion 
of a pressure floor made little difference to the total number or dis- 
tribution of stars formed. 

We acknowledge that as star formation in our simulations may 
happen below the resolution limit, artificial fragmentation can oc- 
cur prior to star formation. These simulations are, therefore, not 
suitable for studying such small scale phenomena as the spatial dis- 
tribution of stars within star forming regions or the mass function 
of molecular clouds. Instead, we focus on the total star formation 
and distribution of star formation above our resolution limit. 



2.6 Cosmological Simulations 

We used cosmological simulations of a dwarf galaxy simulated 
with and without H2 to study the effect of H2 on the star for- 
mation and structure of the stellar disk. Cosmological simulations 
of galaxies are necessary to study H2 throughout galaxy evolu- 
tion and in galaxies with realistically built structure. Dwarf galaxies 
were chosen as an extreme environment for studying gas dynamics, 
star formation, and the effect of SN feedback. Here we describe the 
cosmological simulations presented in this paper. 

We simulated the formation a field dwarf galaxy with an ap- 
proximate final virial mass of 4 X 10"Mq with H2 (DH2) and 
compared it to the same galaxy simulated without H2 (DnoH2). 
We also included two lower resolution versions of the dwarf galaxy 
with H2 (DH2_mr and DH2_lr). While we have verified that our 
H2 formation recipe holds across the range of resolutions of these 
simulations, star formation and feedback are also sensitive to reso- 
lution and dwarf galaxies may be particularly sensitive to changes 
in feedback. We, therefore, included the lower resolution simula- 
tions to ensure that star formation converges at this resolution for 
this mass of galaxy when H2 is included. Aside from differences in 
resolution and as to whether or not H2 was included, all the simula- 
tions are identical. The simulations include cooling from H and He 
ions and from metal lines, cosmic UV background radiation, and a 
blastwave model for SN feedback, as described in section i]2.1l 

The initial c onditi ons for this dwarf galaxy were first used in 
iGovemato et al.l ( I2OICII) as dwarf galaxy 1 (DGl). In that paper, it 
was simulated with neither full metal line cooling nor H2 and pro- 
duced a bulgeless dwarf galaxy. The initial conditions consisted of a 
25 Mpc box surrounding a halo selected from a low-resolution, DM 
only simulation. The initial power spectrum used for the linear den- 
sity field was calculated using CMBFAST code. We assumed a con- 
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Table 1. The cosmological simulations of the dwaif galaxy, listed with their resolution and star formation parameters. The two simulations listed 
in bold are the high-resolution rans and used in the plots throughout the remainder of the paper. The average smoothing length of the gas particles, 
(h), is shown here for the disk gas, which was defined as having a density of 1 amu/cc or greater. All simulations use the same initial conditions 
and include cooling from H and He ions and from metal lines, a cosmic UV background radiation, and a blastwave model for SN feedback. 



Particle Mass [Mc. 



Star Formation Parameters 



Name 



Ha 



Dark 



Star- 



Gas 



Softening 
Length [pc] 



{h) of Disk 
Gas [pc] 



Pmin Tmax 

[amu/cc] [K] 



DnoH2 no 

DH2 yes 

DH2_rm' yes 

DH2Jr yes 



16,000 
16,000 
37,000 
126,000 



1000 3300 

1000 3300 

2400 5900 

8000 20000 



87 
87 
115 
173 



60 
60 
73 
100 



0.1 

0.1 

O.IXh^ 
0.1 Xh, 



100 
0.1 
0.1 
0.1 



2e4 
le3 
le3 
le3 



cordance, flat, A-dq minated cosmology with values from WIVIAP3 
dSpergel et al.l2003h . 

The two high-resolution simulations (DH2 and DnoH2) dif- 
fer both in terms of whether H2 was included and what star for- 
mation recipe was used. In order to study the effect of limiting 
star formation to regions with H2, we 1) multiplied the star for- 
mation efficiency, c*, by the fraction of H2 and 2) lowered both 
the threshold density and temperature, as previously described in 
§[23] Aside from these changes, o ur parameters for star fo rmation 
and fee dback were the same as in lGovemato et al.l (l2010h . We as- 
sumed a lKrouDaetalJ ( ll993l) IIVIF and calculated SN energy using 
the scaling factor Esn ~ 10^^ ergs per SN. The star formation pa- 
rameters and information on resolution for each of the simulations 
can be found in Table[T] 



3 RESULTS 

We analyzed the effects of including H2 and H2 -based star forma- 
tion by comparing the total amount, distribution, and history of the 
gas and stars in our simulations. In order to compare our simula- 
tions to observed galaxies, we generated mock-observations of our 
galaxies in multiple bands using SUNRISE I Jonsson 2006). SUN- 
RISE is a ray-tracing radiative transfer program that includes dust 
scattering and absorption. We used it to measure photometric col- 
ors, luminosities, and star formation indicators such as Ha, FUV, 
and 24 ^m. In addition to simulating the stellar emission from the 
galaxy, we generated HI velocity cubes. These velocity cubes were 
used for calculating the line profile when determining the galac- 
tic rotation speed and for determining the HI surface density when 
calculating the resolved and global K-S relations. We produced the 
velocity cubes by using the smoothing kernel to distribute the gas 
particles spatially. We then calculated the expected amount of 2 1cm 
emission from HI and binned it in x, y, and velocity space. In or- 
der to make these results comparable to observations, we generated 
the velocity c ube a t a resolution and sensitivity typical of THINGS 
JWalter et al.ll2008l) for a dwarf galaxy at a distance of 5 IVIpc. As 
was done for the dwarf galaxies in THINGS, we generated data in 
128 velocity channels with 1.3 km/s bins and a 1.5" pixel size. We 
then spatially smoothed the data cubes using a Gaussian beam with 
a full width at half max of 10"xlO" Finally, we made a sensitivity 
cut and discarded all emission from cells below a 2a noise limit, in 
which a = 0.65 milli-Jansky/beam. 

The final global properties for the galaxies simulated with H2 
(DH2, DH2.mr, and DH2Jr) and without H2 (DnoH2) are listed 
in Table [2] These properties include the virial mass, stellar mass, 
mass of gas within Rvir, and mass of disk gas. The disk gas for this 



puipose was defined as being the mass of HI plus the mass of H2 (if 
applicable) multiplied by a factor of 1 .4 to represent the presence of 
cold He. The values of the Johnson B-band magnitude, r25, and B- 
V color were calculated from SUNRISE observations. The average 
and standard deviation of the SFRs were calculated over the last 2 
Gyrs of the simulation. Finally, the average metallicity of the gas 
(12-1- log(0/H)) was calculated for the star forming gas, as further 
described in ii l3.2l 



3.1 Resolution 

From Table 2, it is apparent that the medium resolution simula- 
tion, DH2_mr, was able to reproduce the star formation of DH2, in- 
cluding the total stellar mass, the stellar distribution (r25), and the 
current SFR (B-V color and (SFR)2Myr)- DH2_lr was a less close 
match to DH2 and produced a slightly brighter and larger galaxy 
with approximately 1.16 times the total mass of stars. We, there- 
fore, conclude that star formation in this galaxy converges for reso- 
lutions of DH2_mr and greater. The initial masses of the gas parti- 
cles in DH2Jr and DH2.mr are 2 x lO" and 5900 IVI©, respec- 
tively. The convergence of star formation between these two mass 
resolutions is consistent with the recommendation from the isolated 
Milky Way-like simulations for gas particle masses of 1.4 x 10^ or 
less (§|2!4j. The total gas mass in the galaxies is more sensitive to 
resolution and there is an approximately 12% difference in the gas 
mass within the virial radius between DH2_mr and DH2_hr. This 
change indicates a slight decrease in the amount of gas within the 
halo with increasing resolution. 



3.2 Stellar and Baryonic IVIass 

In order to compare the global properties of the simulations, we 
first studied the total amount of stars and gas contained in a DM 
halo. In dwarf galaxies, being able to correctly simulate the lumi- 
nosities of galaxies as a func tion of their DM h alo mass is key to 
the missing satellite debate ( lMooreetal1ll999t) . This debate ad- 
dresses why simple extrapolations of the luminosity function from 
DM-only simulations predict more low-mass satellites than are ob- 
served. The luminosity of the galaxy and the amount of baryons 
contained in a DM halo are a function of the efficiency of star for- 
mation, the amount of gas that is accreted, and the amount of gas 
that is lost as the result of feedback. The inclusion of H2 has the 
potential to affect the amount of gas that turns into stars and the 
ainount of SN feedback. The efficiency of SN feedback is affected 
by the state of gas, with clumpier gas leading to greater gas loss. 
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Table 2. Final properties of tlie dwaif galaxy simulated with and without H2. The two simulations listed in bold are the high resolution simulations 
and are used in the plots throughout the remainder of the paper. Disk gas mass was defined as 1.4(HI + H2), as is typically done in observations to 
account for the mass of cold He. The average and standard deviation of the SFR were calculated over the last 2 Myr of the simulations, using bins 
of 5 X lO'^ years. The metallicity, 12 + log(0/H), was calculated for the star forming gas, as described in 13. 31 



Mass [IO^Mq] 



Name 


Virial 


Stellar 


Gas 


Gas 


Mb 


^25 


B-V 


(SFR> 


0"SFR 


12 -1- log(0/H) 








(/? < Rvir) 


(Disk) 




[kpc] 




[M0j/r-i] 






DnoH2 


360 


2.1 


17 


4.3 


-15.27 
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29 


8.0 
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0.26 
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Therefore, the inclusion of H2 could change the amount of gas lost 
by producing a denser, colder ISM. 

In order to measure the effect on the stellar and baryonic mass, 
we determined the dark, stellar, and gas mass of the galaxies (Ta- 
ble|2]l and calculated their positions on both the standard and bary- 
onic TuUy-Fisher relationships. While the standard Tully-Fisher re- 
lationship relates the luminosity of galaxies to their rotation ve- 
locity, the baryonic Tully-Fisher relationship relates the amount of 
baryonic mass in a halo to its rotation velocity. The baryonic Tully- 
Fisher relationship is a particularly important relationship for low- 
mass galaxies because of the wide spread in their stellar-to-mass 
ratios. It has been observed to be an exceptionally tight relationship 
( lMcGauglj[2003) and can, therefore, provide a stronger constraint 
than the standard Tully-Fisher relationship. 

In Figure |6l we compare the luminosities and baryonic frac- 
tions ofour_gal^xi£s_to the observational data compiled in Figure 
6 of lGeha et all (|200^. For both relationsh ips, we measured the 
HI line width at the 50% level (W21) as in lHavne s et al. 1 1999), 
using the 21cm line from a simulated HI velocity cube with the 
galaxy oriented at 45°. When comparing our galaxies to the stan- 
dard Tully-Fisher relationship, the I-band magnitude was calcu- 
lated from a SUNRISE simulated observation. When comparing 
the galaxies to the b aryonic Tully-Fish er relationship, we followed 
the same method as iGeha et alj JioO^ to find the baryonic mass: 
the total gas mass was assumed to be 1 .4 times the HI mass and the 
stellar mass was determined from the SDSS /-ban d magnitudes and 
g - r colors (calculated with SUNRISE) and the lSell et all j2003l) 
mass-to-light ratios. Both data points are within the observational 
spread, although the presence of H2 resulted in a galaxy richer in 
both stars and cold gas. 

As seen by the positions of the galaxies in the Baryonic Tully- 
Fisher and quantified in Table|2l the inclusion of H2 increased the 
baryonic mass to 1.28 times the baryonic mass in DnoH2. Most of 
this increase was in the form of gas: the mass within Rvir for DH2 
was 1.29 times higher than in DH2. The stellar mass increased to 
only 1.16 times higher than in DnoH2, indicating a slightly lower 
stellar mass to gas mass fraction. 

We investigated the source of DH2's larger baryonic fraction 
by finding the total mass of gas ever in the galaxies and ever ejected 
by SNe from the galaxies. We calculated the mass of gas that was 
ever in the galaxy by identifying the all gas particles that were in- 
cluded within the main halo at any timestep. We then found the 
mass of gas lost from the galaxy by identifying all gas particles 
that had ever been heated by SNe and that were once part of the 
main halo but were outside the halo at z = 0. Our method for cal- 
culating the mass of gas ever in the galaxy does not include gas 
that became stars in separate smaller halos prior to their accretion. 
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Figure 6. Standard (left) and baryonic Tully-Fisher ( right) relations for the 
two simulations, overlaid on observational data from lGeha et alj j2Q06l) for 
dwarf galaxies (grey filled circles). The diamond represents DnoH2 and the 
square represents DH2. Both simulations lie within the observed scatter at a 
redshift of zero but the inclusion of H2 resulted in a slightly higher stellar 
and baryonic mass. 



It also has the potential to miss gas that was accreted onto and ex- 
pelled from the halo between two adjacent timesteps. However, we 
are interested in the comparison between the two simulations rather 
than the absolute amounts so these systematic uncertainties are rel- 
atively unimportant. 

We found little difference in the total mass of gas accreted onto 
DnoH2 and DH2. DnoH2 and DH2 had a total masses of gas ever 
in the galaxy of 3.71 x 10^ M© and 3.64 x lO' M©, respectively 
There was actually slight increase in the amount of gas expelled 
by SNe beyond the virial radius in DH2 compared to DnoH2. In 
DH2, a total of 5.4 x 10* Mq was expelled by SNe, compared to 
3.9 X 10** M0 expelled from DnoH2. The decreased baryonic mass 
of DnoH2 compared to DH2 was, therefore, not the result of either 
lower amounts of gas accretion or greater amounts of SNe-driven 
gas loss. Instead, the greater amounts of gas loss in DnoH2 came 
from gas being smoothly stripped from the outskirts of the halo. 
This stripping of gas may be an indication that gas in DnoH2 was 
less tightly bound than in DH2 or it may be the result of stochastic 
variation in the mergers. 
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3.3 Metallicity 

The total stellar mass, in addition to the amount of gas lost and ac- 
creted, is responsible for the metallicity of the galaxy. Given the 
importance of metallicity to the fomation and shielding of H2, 
we show the simulated galaxies on the observed mass-metallicity 
relationship ( Figure [7}. In thi s figure, the observational data is 
taken from the Lee et al. 72006*) set of dwarf irregular galaxies and 
the [Xremonti et al. (2004) fits for their survey of SDSS galaxies. 
The latter is scaled down by 0.26 dex to reflect the tendency of 
the method use d to overestimate the oxygen abundance relative to 
other methods jTremonti etal1l2004 lErb et allliooel brooks et al] 
I2OO7I) . DH2 lies o n top of an observed galaxy and very close to 
the Tremonti et al.l 12004) fit. The metallicity of DnoH2 is slightly 
lower than that of observed galaxies of the same stellar mass but 
may lie within the spread predicted from lTremonti et al1(l2004h . 

The stellar masses of the simulated galaxies used in the mass- 
metallicity relationship were calculated from mo ck-observa t ions o f 
the stellar emission using the same scaling as in iLee et al.l ( |2003) . 
The mass-to-light ratio was determined from the B-K color and 
combined with the magnitude in the K-band to produce a stel- 
lar mass. The metallicities in this plot were determined for the 
star forming gas. In order to account for the fact that observations 
of gas metallicity are determined for HII regions, the metallici- 
ties of the g as we re sc aled by the star formatio n probability, as in 
iDave et alj ilOOH) and 'Finlato r & Davd ( |2008|) . Compared to the 
average metallicity of all HI and H2 gas within the radius where 
star formation has occurred in the last 100 Myrs, this scaling has 
the potential to change the metallicity by up to 0.1 dex. In both 
these galaxies a substantial amount of low-metallicity HI gas lies in 
an extended disk just beyond the star formation radius. Therefore, 
while the star formation takes place in primarily higher metallicity 
gas, extremely low-metallicity gas chemistry remains important to 
the galaxy at z = 0. 



3.4 Star Formation Histories 

In order to determine what circumstances led to differences in the 
stellar masses of the two simulations, we compared the evolution 
of the star formation his tories of DH2 and DnoH2. Previous au- 
thors (e.g. R obertson & Kravtsovll20o3: iGnedin & Kravtsovll2010l : 
iKrumholz & Dekelll201lh have suggested that introducing a metal- 
licity dependency to star formation by linking star formation to H2 
could suppress star formation at high redshift. Using galaxy evolu- 
tion si mula tions with H2 -based star formation. iGnedin & Kravtsovl 
( l2010h and lKuhlen eTai] |201 1) hove show reduced star formation 
up to a redshift of 3 and 4, respectively. 

We examine whether low-metallicities at high redshift may 
have suppressed star formation in DH2 relative to DnoH2 by com- 
paring the evolution of the gas-phase metallicity to the SFR and 
total stellar mass (Figure |8). The top panel shows the mean gas- 
phase metallicity of the star forming gas in the largest progenitor 
as a function of time. Notably, there is almost no evolution in the 
gas-phase metallicity since a redshift of four. For the majority of 
the galaxies's histories, the inflow of metal-poor gas has balanced 
the production of metals from stars. The second and third panels 
show the cumulative star formation histories and the star formation 
histories of the galaxies. Unlike the metallicity panel, these panels 
show the history of all star particles that are in the galaxy at 2: = 0, 
rather than only those particles that were in the galaxy at a given 
redshift. 

Figure [8] shows slight evidence for a decrease in high redshift 
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Figure 7. Mass-metallicity relationship with DnoH2 (diamond) and DH2 
(square) overlaid on observational data. The stellar masses of the simu- 
lated galaxies were calculated in an observationally-motivated fashion from 
the Johnson-Morgan K-band luminosity. The metallicities shown for the 
simulated galaxies ai'e the mean metallicities of the gas particles scaled by 
the probability of star fonnation . The filled grey circles represent observed 
dwarf irregular galaxies from Lee et al. (2006). The grey lines are the obser- 
vational fits from Tremonti et al. (2004): the solid line is the best fit line, the 
dot-dashed line is that fit extended to lower stellar masses, and the dashed 
and dotted lines show the 1 and 2 cr spread, respectively. The metallicity 
of DH2 is consistent with the observed galaxies. The metallicity of DnoH2 
is slightly lower than the observed g alaxies but may be co nsistent with the 
spread in metallicities determined in lXremonti et al.l | |2004|) . 



star formation. DH2 does initially have a slightly lower stellar mass 
than DnoH2 but hy z — 4, DH2 has begun to surpass the stellar 
mass of the other galaxy. At this same redshift, the metallicity of 
the star forming gas settles to a near constant state. With our lim- 
ited number of high redshift outputs, it is not clear if the decreased 
stellar mass at z > 4 is the result of a higher star formation thresh- 
old or simply stochastic variation in the merger history. Regardless, 
any suppression of star formation compared to DnoH2 is of brief 
duration. 

In general, star formation in our simulations appears very ro- 
bust to the specifics of the star formation recipe. Essentially, if the 
star formation efficiency is decreased either artificially or because 
of decreased H2 content, the star formation may be delayed but the 
particle will continue to grow denser and more molecular and will 
eventually form a star. For these simulations, SN feedback is the 
primary regulator of star formation. The relative importance of H2 
formation and SN feedback and how it compares across different 
simulations is addressed in further detail in the discussion. 



3.5 Color and Extent of Stellar Disk 

In addition to changes in the total mass, we examined the effect 
of H2 on the current star formation. From the data in Table |2] the 
SFR over the last two Myrs was four times higher in DH2 than in 
DnoH2. This increase is apparent from an observational standpoint 
in the colors of the galaxies. In Figure |9] we compared the B-band 
magnitudes and B-V c olors of the gal axies to observational data 
for dwarf galaxies from Ivan Zeel(l2000l) . While both the colors and 
magnitudes of DH2 and DnoH2 are consistent with the observed 
data, we found that the addition of H2 resulted in brighter and bluer 
galaxies. The brighter B-band magnitudes of DH2 were the result 
of both a larger stellar mass (ij |3.2t and more recent star formation. 
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Figure 8. The evolutions of the gas metaUicity and stellar mass, and the star 
formation histories. In each of the panels, the dashed line represents DnoH2 
and the solid line represents DH2. The top panel shows the mean gas-phase 
metallicity for the largest progenitor of final galaxy. The middle and bottom 
panels show the star formation histories. In these two plots, all star particles 
within the galaxy at z = are included. The middle panel shows the cumu- 
lative star formation histories (stellar mass formed as a function of time). 
The bottom panel shows the SFRs over time. There is a very slight decrease 
in star formation in DH2 prior to z = 4. In general, though, the star fonna- 
tion rates track together with the slightly higher rates in DH2 reflecting the 
greater disk mass. 
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Figure 9. B-V color versus B-band magnitud e for the simula ted galaxies 
overlaid on observations of field dwarf's from Ivan ZeS fcOOOl) (grey aster- 
isk). The diamond represents DnoH2 and the square represents DH2. The 
inclusion of H2 resulted in the dwarf galaxy being brighter and bluer, indi- 
cating a greater amount of on-going star formation. 
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Figure 10. Color as a function of radius for DnoH2 (dashed line) and DH2 
(solid line). The vertical lines represent r25 for each galaxy. The bluer col- 
ors at the optical radii in DH2 indicate on-going star fonnation throughout 
a larger area of the disk. 



with greater amounts of star formation at z = 0, but also caused star 
formation to occur over a larger extent. 



In order to compare the distribution of stars and star forma- 
tion in the two galaxies at z = 0, we used SUNRISE to determine 
the radially-binned B-V colors for the galaxies when oriented face- 
on (Figure [Tot. Overlaid on the plot are vertical lines representing 
the locations of r25. DH2 had much more spatially extended star 
formation, as can seen by the blue colors that extend out to sev- 
eral kpc in radii. This extended star formation is evident in a larger 
stellar disk, characterized by the greater value of r25. This distribu- 
tion of colors is consistent with obse r vation s of the color profiles of 
dwarf irregular galaxies in Ivan ZeeldlOOOl) . In contrast, DnoH2 is 
distinctly blue at the center and surrounded by a dim, redder pop- 
ulation. This color profile indicates that the small amount of star 
formation at z = was extremely concentrated at the center of the 
galaxy. Therefore, the addition of H2 resulted not only in galaxies 



3.6 Kennicutt-Schmidt Relation 

Based on the changes to the star formation caused by the inclu- 
sion of H2, it is important to compare the star formation laws in 
both simulated galaxies to the observed K-S relation. The K-S re- 
lation has be en observed both as the global K-S relation on galaxy 
wide scales l lKennicut3l 19981) and as the resolved K-S relation on 
hundred-parsec scales (^Big iel et al .I2OO8I) . The existence of a global 
K-S relation implies a consistent relationship between star forma- 
tion and gas surface density when averaged over the entire disk of 
the galaxy. On the sub-kpc scale, the correlation requires this re- 
lationship hold even as smaller star-forming regions begin to be 
resolved. 

Star formation recipes in simulations that are based on the 
free-fall time, as ours is, imply a K-S relation for gas above the star 
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formation density thresiiold (pmin)- Such recipes tiiat allow star 
foiTnation in all disk gas by having a low pmin, such as 0.1 amu/cc, 
are more directly related to the global K-S relation. In contrast, 
star formation recipes, such as ours, which limit star formation to 
higher density regions through the use of high-density thresholds or 
H2-based star formation are more directly related to the resolved 
K-S relation. Resolving high-density regions in simulations and ty- 
ing star formation to those regions through the use of a high den- 
sity th reshold improves the spatial distribution of stars dSaitoh et al.l 
1200^ ). The use of a high-density threshold, however, removes the 
direct connection between the SFR and the global gas distribution. 
Fitting the global K-S relation requires that the density structure 
of the gas result in the appropriate fraction of high density, star- 
foiTning gas. In order to verify that star fomation in our galaxies 
is occurring at the coiTect rate over all spatial scales, we replicated 
both the resolved K-S relation and the global K-S relation. 

In order to replicate the resolved K-S relation, we determined 
HI, H2, and SFR surface densities in 750 pc x 750 pc boxes within 
r25 and compared them to the gas and SFR surface densities mea- 
sured over the sam e spatial scale and within r25 for dwarf galaxies 
( iBigiel et alj201(]|) . Our HI surface densities were determined from 
the zeroth moment of a simulated velocity cube of 21cm emission 
with average THINGS resolution for the galaxy oriented at a 45° 
angle. This velocity cube is described in detail at the beginning of 
§3. We corrected for the inclination by multipl ying all surface den - 
sities by sin(45° ). We calculated the SFRs as in lBigiel et aD j2008l) , 
from a combination of FUV and 24/im flux in a SUNRISE mock- 
observation. The surface densities of H2 were determined directly 
from the simulation, without reference to CO. 

The observation al K-S data we co mpared to is a subset of 
th e data presented inlBigiel et al.l ( I2OI0I) . The K-S data presented 
in lBigieletal.1 feOlOh is for all the THINGS galaxies, including 
the more metal-rich spiral galaxies and the low-metallicity dwarf 
galaxies. Since the surface density at which HI transitions to H2 
is metallicity dependent, we limited our comparison to only the 
metal-poor dwarf galaxies from THINGS: DDO 53, DDO 153, 
Holmberg I, Holmberg II, IC 2574, MSI Dwarf A, and M81 Dwarf 
B. Thes e galaxies range i n metallicity from 7.54 <I2 -1- [O/H] 
< 7.94 l lMoustakajfioci^) and the majority of points come from 
IC 2574, which has a metallicity of 12 -1- [O/H] = 7.94, very similar 
to DH2. The gas surface densities shown are equal to 1.4Ehi; H2 
is not included in the observations because of the lack of observed 
CO and uncertainty in Xco. Low metallicity dwarf galaxies are HI 
dominated, however, so ne glecting H2 should have only a small 
effect on the surface density. Isolatto et al .1 ( 120 II I) (not shown in the 
figure) is another recent observational study of the resolved K-S 
relation in low metallicity environments. In that study, very high 
resolution ( 12 pc) measurements of the HI and total gas surface 
density were taken in order to calculate the H2 surface density in- 
dependent of CO emission. They observed similar v alues of Ssfr 
at 3 times higher gas surface densities than the Bigie l et al] j20IOl) . 
This discrepancy may be due to physical differences between indi- 
vidual galaxies or the observational techniques used. 

Both simulations reproduce the steep slope of the resolved K- 
S relation at low surface densities. In DnoH2, the steep slope at 
lower surface densities was generated by the minimum density re- 
quirement for star formation. In DH2, the steep slope marks the 
transition from atomic to molecular gas. DH2 has a higher SFE 
than DnoH2, as evident by the lower surface densities at which star 
formation takes place. The higher SFE in DH2 could be because 
the likelihood of star formation increased gradually with increas- 
ing H2 fraction, rather than turning on sharply at the star formation 
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Figure 11. The resolved K-S relation within r25 for the simulations at z = 
0. The simulation data are ov erlaid on obs ervational data of low metallic- 
ity dw arf galaxies from Bigiel et aP j20ld !) and of the SMC lBolatto et al] 



ity dw i 

( 2011). The surface density of HI was computed from a mock-THINGS ob- 
servation and the SFRs were deteiTnined from SUNRISE-generated mock 
observations of the FUV and 24/im flux. The diamonds are for DnoH2 
while the squares are for DH2. Both simulations reproduce the steep slope 
at low surface densities. DH2, however, has star formation at sightly lower 
surface densities, more star-forming gas, and a greater spread of data points. 



threshold density, p — 100 amu. It also could reflect greater het- 
erogeneity in the gas density in DH2 (see discussion), such that the 
high density of the star-forming gas was tempered by low density 
nearby gas. 

In addition to the slightly lower gas surface densities in DH2, 
the primary differences between DnoH2 and DH2 are in the spread 
and number of points. In DH2, the greater amount of spread is 
likely because the density threshold for star formation in DnoH2 
was replaced by the more gradual transition as the gas become 
molecular. The greater number of data points for DH2 is a reflec- 
tion of the greater area in the galaxy forming stars at a redshift of 
zero, as discussed in i]3.5l 

In order to compare DH2 and DnoH2 to the global K-S re- 
lation (Figure II2t , we used SUNRISE-generated Ha emission to 
find the SFR at z=0. We calculated the total HI and H2 surface 
density within r25, where r25 was calculated from a Johnson B- 
band SUNRISE observation and HI and H2 surface densities were 
determined as for the resolved K-S relation. We then compared 
the location of the simulated g alaxies on th e global K-S relation 
to both spiral i Kennicut] 1 1 9981 ) and dwarf ( Rovch owdhurv et al.l 
2009) galaxies. Despite the fact that the final SFR in DH2 was four 
times higher than in DnoH2 (Table |2ll, both galaxies have similar 
locations on the K-S relation. This is because the star formation 
in DH2 was much less spatially concentrated than in DnoH2. As 
in the resolved K-S relation, DH2 fits particularly well with the 
observed galaxies while DnoH2 has a slightly higher gas suiface 
density than observed galaxies with similar star formation surface 



© 0000 RAS, MNRAS 000, 000-000 



14 C. Christensen et al. 



densities. Both simulated galaxies are below the K-S relation best- 
fit line and lie within the scatter of observed galaxies. 



It is gratifying that both DH2 and DnoH2 fit more closely with 
the observed dwarf galaxies than spiral galaxies. The lower ratios 
of EsFR to Egas in dwarf galaxies indicate smaller fractions of 
star-forming gas. These lower amounts of star-forming gas are pre- 
sumed to be the result of lower-metallicities, which raise the sur- 
face density where H2 is found, and smaller fractions of gas tha t 
reach high surface densities. While iRovchowdhurv et al.l ( l2009h 
does not contain metallicity measurements for these galaxies, they 
can be assumed to lie on the stellar mass-metallicity relationship. 
With B-band magnitudes ranging from -14.90 to -10.78, they are 
even dimmer than our simulated galaxies (AfB,DH2 = —15.70 and 
MB,DnoH2 = —15.27) and likely have metallicities that are sim- 
ilar or even smaller. Both DH2 and DnoH2 have lower SFEs than 
would be predicted from the global K-S relation because 1) their 
high star formation thresholds limit the fraction of gas capable of 
forming stars and 2) these galaxies have fractions of star-forming 
gas more similar to observed dwarf galaxies than spiral galaxies. 
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In addition to calculating the location of the galaxies on the KS 
law, we calculated the molecular gas depletion time (E H2 /Esfr) 
on 750 pc scales. H2 forms at the correct surface densities and 
the star formation efficiency has been tuned to produce appropri- 
ate amounts of star formation. However, the inability of the code 
to resolve the interior structure of molecular clouds in connection 
with stellar feedback results in artificially low amounts of H2 near 
recent star formation. Essentially, in our model it is impossible for 
H2 to remain in a gas particle that has recently formed a star be- 
cause of the spatial scale of SN feedback. We account for this effect 
when calculating E H2 by adding in the amount of H2 when star 
formation occurred. Using this rough estimation, we find a smaller 
average molecular depletion time of 10* years, compared to molec- 
ular depletion time on 100 parsec scales of 1.6 x lO" years for the 
SMC ( Bolatt o et al, jOl 1) and 2.0 x 10^ years for the THINGS spi- 
ral galaxies ( iBigiel et alj|2008h . Note that the observed molecular 
depletion times are generally calculated for higher star formation 
rate surface densities: 10~^ < Esfr < 10^^, as opposed to our 
data which had Esfr < 10"'^. The shorter depletion times of our 
simulations may stem from the uncertainty in our calculation and 
the different Esfr of the observations it is being compared to or 
it may indicate a need to raise the clumping factor, Cp, used in the 
H2 formation calculation. 



The lack of contrast between DnoH2 and DH2 in either ver- 
sion of the K-S relation is notable. Even with very different mod- 
els for the interstellar media and for star formation, the results for 
both are similar to each other and to the observed data. For these 
two simulations, the relationship between gas surface density and 
star formation is similar since both star formation recipes are based 
on the free-fall time and both impose a high-density threshold (in 
DnoH2, an explicit threshold, in DH2, an implicit one based on the 
H2 abundance). The similar results for DnoH2 and DH2 indicates 
that while the K-S relation is important to reproduce, it is an im- 
precise tool for distinguishing between these two simulations. A 
similar correspondence between the K-S relation generated with a 
high density threshold for star fo rmation and a H2 dependency was 
presented in iKuhlenetalJ j201lh . 



Figure 12. The global K-S relation for the simulati ons at z = ove rlaid 
on observational data from spiral galaxies (asterisks; * Kennicut3ll99 8') and 
dwarf galaxies (filled circles; Roychowdhury et al. 200^^)~Thediamond rep- 
resents DnoH2 and the squai'e represents DH2. The cross at the top of the 
plot represents the average en'or for the observational data in (Kennicut] 
[1998). The line is the K-S power-law with an exponent of 1 .4 for spiral and 
starbursting galaxies. DH2 fits well with the observed dwarf galaxies and 
DnoH2 lies very neai' to the observed galaxies. 



4 DISCUSSION 

4.1 Changes in Gas Thermodynamics 

The increase in baryonic mass and star formation at z = is most 
likely the result of the increased amount of cold, dense gas formed 
in our simulation, caused by the inclusion of H2 and the shielding 
and cooling that accompany it. In incorporating H2, we included 
shielding of HI by dust, in addition to dust and self-shielding of 
H2. This shielding of HI and H2 gas from photoheating drastically 
reduces the heating rates of the dense disk gas. In simulations that 
have both shielding and low temperature (T < 10 ) cooling, either 
from metal lines of H2 (see the following paragraphs), the lack of 
heating will result in cold gas. In DH2, both of these conditions 
are met, whereas in DnoH2, the gas was unshielded, in addition to 
having smaller cooling rates. 

Evidence for the increased amount of cold, dense gas caused 
by our H2 implementation is apparent when comparing the star- 
forming gas in DnoH2 and DH2. Figure[T3]contains histograms of 
the density and temperature of the gas from which stars formed in 
the simulations. In simulations with H2, the formation of a cold 
phase of the ISM is apparent in the colder and denser star-forming 
gas. The star-forming gas in these simulations is much more similar 
to the star-forming gas in actual galaxies. 

While the incorporation of shielding is primarily responsible 
for the increased amount of cold gas, H2 also acts as an additional 
source of cooling at temperatures between 200 and 5000 K (Fig- 
ure [T4j- The importance of cooling from H2 in relation to other 
sources of cooling was heightened in the dwaif galaxy in part be- 
cause of its low metallicity. The metallicity of the star forming gas 
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Figure 13. Normalized histogram of the densities (left) and temperatures 
(right) of the gas particles that formed stars. The dashed line represents 
DnoH2, while the solid lines are for simulations with H2 at the three dif- 
ferent resolutions, DH2Jr (thin), DH2.mr (medium weight), DH2 (thick). 
The densities and temperatures of gas that form stars in simulations with 
H2 are very similar across all resolutions. Star fomiation in simulations 
with H2 in higher density and colder gas than in DnoH2. 
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Figure 14. Cooling rate for gas at a redshift of zero. The left panel com- 
pares the rates of gas cooling between DnoH2 (black) and DH2 (grey). The 
points represent individual SPH gas particles. The metallicity of gas parti- 
cles shown in this image covers the range —1.5 < log(Z/ZQ) < —0.5. 
When H2 was included, the gas was able to cool to lower temperatures 
and had a higher cooling rate at temperatures below 10* K. The light panel 
compares the rates of gas cooling in DH2 (grey) to the estimated maximum 
cooling rate from CII, assuming all carbon is signally ionized. The esti- 
mated CII cooling rate is similar to that from H2 in low-metallicity regions 
with molecular gas. 



in DH2 was 12 -1- log(0/H) = 7.9 a.nd the metallicity of the disk 
gas beyond the star forming radius was even lower. At these metal- 
licities, the metal line cooling will be less effective than in spiral 
galaxies. 

The relative importance of H2 to metal line cooling in our 
simulations is also heightened by the different physics assumed 
when calculating the H2 and the metal line cooling, leading to 
an underestimate of the metal line cooling rates compared to the 
H2 cooling rates. The calculation for the H2 abundance includes 
a model for shielding of HI and H2 and a model for dissociat- 
ing radiation from the interstellar radiation field (ISRF). The metal 
line cooling rates were introduced to GASOLINE in order to bet- 
ter model the enr ichment and cooling of the intergalactic medium 
dShen et alj I2OIO ). These cooling rates were calculated assum- 
ing optically thin gas for a range of temperatures, densities, and 
redshift-de pendent UV back ground radiation using CLOUDY (ver- 
sion 07.02: lFerland et al.ll998il . These assumptions are valid in the 
intergalactic medium. However, in the ISM, the lack of UV radi- 
ation from the ISRF, even with the inclusion of cosmological UV 
background radiation, results in a much lower UV flux. The lack 
of UV radiation decreases the rate of cooling from forbidden line 
transitions in unshielded gas, in particular CII. With more forbid- 
den line transitions, the cooling rates for the low temperature gas in 
DnoH2 would be expected to substantially increase. The observed 
difference in cooling rates between DH2 and DnoH2 for gas with 
temperatures less than lO^K is, therefore, exacerbated by the un- 
derestimate of metal line cooling. 

We estimated the maximum amount of cooling that could the- 
oretically be produced by CII and compared it to the cooling rates 
of the gas in DH2 (Figure [74t. To make this estimate, we assumed 
that all carbon in the gas below 10*K was in the CII state and 
that CII cools at a r ate of 2.89 x 10^^" N(CII)/N(H) erg s'^ H'^ 
dLehner et al]|2004h . We calculated the amount of carbon from the 
oxygen abundances for ea ch gas particle using the equation for C/0 
from .Gamett et al] d 19951) . At the low metallicities of the gas, we 
found that the estimated maximum amount of cooling from CII is 
approximately equal to the cooling from H2 in regions with molec- 
ular gas. 

In the future, extending our calculations of the ISRF and 
shielding to include chemicals other than H2 will allow us to better 
model the cooling of dense gas in the ISM. At present, though, the 
introduction of H2 in combination with shielding is an important 
step in the simulation of the cold ISM. The addition of shielding 
and efficient low temperature cooling will quickly lead to the cre- 
ation of cold gas in regions shielded from photoheating, regardless 
of the source or exact rate of that cooling. Essentially, so long as 
heating is reduced by shielding, the presence of low-temperature 
cooling will lead to the creation of cold gas. Therefore, to the ex- 
tent that the cold gas is resolved, this implementation of H2 leads 
to a good first approximation of what the ISM would be like in sim- 
ulations with more advanced models of low-temperature cooling. 

The distribution of dense gas in both DH2 and DnoH2 can 
give insight in to how changes to the thermodynamics affected the 
amount and distribution of star formation. In DnoH2, star forma- 
tion was mostly confined to the center of the galaxy by a redshift of 
zero. In contrast, when H2 was included star formation continued 
to occur throughout the disk. As can be seen in Figure[T5] this effect 
does not seem to be due to higher average surface densities at large 
radii. The surface density of DnoH2 is, in fact, slightly higher than 
in DH2 within r25 . However, individual gas particles in DH2 reach 
high densities at larger radii than in DnoH2. The spikes in FigurefTSi 
represent clumps of dense, cold gas in the disk. This clumpier ISM 
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Figure 15. Gas density (left vertical axes and points) and radially averaged 
surface density (right vertical axes and lines) as a function of radius for 
DnoH2 (top) and DH2 (bottom). The radially averaged surface density is 
similar in both galaxies and slightly higher in DnoH2. The spread in densi- 
ties at a given radii is higher in DH2 than in DnoH2, indicating a clumpier 
ISM. 

was the result of the formation of cold gas in DH2. In that simu- 
lation, high-density clumps were able to form throughout the disk, 
whereas in DnoH2 the only area with high density gas particles was 
at the center of the galaxy. 

The clumpier ISM in DH2 also affected the amount of star 
formation at z = 0. Because more of the gas in DH2 was capable 
of star formation than in DnoH2, star formation continued to occur 
at a higher rate. Broad conclusions about the effect of the inclusion 
of H2 on star formation in galaxies of different masses and envi- 
ronments cannot be drawn from one simulation. However, this test 
galaxy illustrates how the incorporation of H2 results in a clumpier 
ISM and the level of importance that effect can have to modeling 
star formation in general. In the future, we will increase our set 
of simulated galaxies in order to study H2 within the context of 
different galactic properties. 

In addition to changing the distribution of star formation at z 
= 0, the ability of gas to reach colder temperatures in DH2 is also 
potentially responsible for the increased mass in the disk of the 
galaxy. Both DnoH2 and DH2 accreted the same fraction of their 
gas over their lifetime and a slightly higher fraction of gas was ex- 
pelle d from DH2. This increased gas loss is the effect of a clumpier 
ISM taovemato et"ai]|201(]| : iBrook et al.ll201ol I2OIII) . However, a 
greater mass of gas was stripped form the halo of DnoH2 than in 
DH2, possibly because the greater ability of gas in DH2 to cool led 
to the gas being more tightly bound. 

By linking star formation of H2, a metallicity dependency is 
introduced into the star formation law. As metallicity decreases, the 
lack of shielding causes H2 and, consequentially, star formation to 
take place in higher density gas. The result is an effective star for- 
mation threshold that shifts with metallicity - someth ing that fixed 
density thresholds cannot produce dKuhlen et alj201lh . We observe 
the Hl-to- H2 transition to shift to higher density with decreasing 
metallicity in our model of an isolated MW galaxy. Similarly, the 
resolved K-S relation for the simulated dwarf galaxies fits with that 
of low-metallicity dwarf galaxies. Finally, we find that the simu- 
lated dwarf galaxies, like observed ones, lie below the global K-S 
relation for disk galaxies, indicating a decreased star formation ef- 
ficiency on global scales compared to disk galaxies. This decrease 
is the result of smaller fractions of star forming gas, either because 



the lower metallicity causes an effective increase in the star forma- 
tion threshold or because these galaxies have a lower fraction of 
high surface density gas. 

We found a lower fraction of the baryonic mass in the form 
of stars in DH2 but the overall increase in the amount of dense, 
disk gas led to an increase in the global SFR. While DH2 had a 
slightly decreased star formation efficiency, the increased amount 
of baryons led to a somewhat higher stellar mass. At high red- 
shift when the metallicities were even lower, DH2 forms only 
slightly fewer stars than DnoH2. This resu lt initially appears c on- 
tradictory to the sim ulations presented in iKuhlen et al.l ( l201lh . In 
iKuhlen et alj j201 ll) . equilibrium H2 abundances were used to cal- 
culate star formation in simulations of dwarf galaxies evolved to 
z = 4, resulting in dramatically decreased stell ar mass. The sim u- 
lations presented here differ from those in Kuhl en et al.l 1 I2OI ih in 
the use of a non-equilibrium H2 calculation, a lack of a minimum 
H2 fraction for star formation, which allows for star formation in 
a greater variety of environments, and the shielding of H2 and HI 
gas from photoheating. We also included a more efficient model 
for SN feedback, in comparison to Kuhlen et al. (2011). The latter 
differe nce is a particularly im portant. 

In lKuhlenetai] j201lh efficient SN feedback was intention- 
ally not included in order to isolate the effects of an H2 -dependent 
threshold. In our simulations, the efficient feedback is included to 
more closely simulate the behavior of gas in actual galaxies. Such 
feedback is responsible for removing baryons from the galaxy, re- 
sulting in an appropriate fraction of baryon mass in the disk and 
an appropriate central concentration of the galaxy. The result of in- 
cluding efficient feedback is that the feedback is by far the strongest 
regulator of star formation. In our simulations, the balance between 
the heating of gas by feedback and the cooling of gas in the galaxy 
is more important to setting the global SFR than the change in 
star formation efficiency produced by an H2 -based star formation 
recipe. We, therefore, most strongly see the effect of H2 on star 
formation indirectly through how it influences the amount of gas 
in the disk and where the cold, dense star-forming gas is found. 
Given the differences in H2 model, star formation, thermodynam- 
ics, and feedback, it is unsurprising that we do not witness the same 
dramatic decrease in star formation. 

Basing our star formation recipe on the abundance of H2 
relies on the premise that the transition to H2 is necessary be- 
fore star formation can begin. The idea that the formation of 
H2 is a necessary precursor to star formation has been driven in 
part by the large number of observations of spiral galaxies link- 
ing CO emission (and by extension, H2) to star formation and 
has been assumed or argued in a number of theoretical papers 
(Elmegreen 2007; Robertson & Kravtsov 2008: Gnedin et a l.l2()09l : 
iPelupessv & Papadopoulos 2009 ; K rumholz et al. 2009b,). How- 
ever, others suggest that H2 is not necessary for star formation but 
instead merely correlated with the conditions that lead to star for- 
mation in normal spiral galaxies. In this model, star formation is 
initiated i ndependent of whether the hydrogen is atomic or molec- 
ular ( e.g. iLi etal.ll200l lOstriker et ailboiol ; iMac Low & Gloveil 
I2OIOI) . Meanwhile, the same processes that trigger star formation, 
such as converging gas flow, cloud collisions, and large-scale grav- 
itational collapse, would generally also result in H2 in spiral galax- 
ies. 

Assuming, though, that star formation occurs because of 
changes in the gas related to the formation o f H2, the que stion 
rem ains as to why this transition is important. Schave (2004) and 
iKrumholz et alj j201 ih argue that the fundamental step that leads 
to star formation is the transition to lower temperatures that hap- 
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pens concurrently with the transition from atomic to molecular hy- 
drogen. In a similar model, Glover & Clark (2011) argued that the 
formation of gas shielded by dust is what leads both to the forma- 
tion of H2 and to the formation of cold gas that leads to star for- 
mation. While recognizing the importance of this debate, we have 
proceeded with our research under the assumption that H2 or the 
shielding and low temperatures that accompany it are intrinsically 
linked to star formation. In the future, the ability to alter the depen- 
dency of star formation on H2 in models of galaxy formation could 
be important for testing the relationship between the two. 

A second area of research that warrants further investigation 
is the connection between H2 and carbon monoxide, CO. Because 
of the faintness of H2 emission at the temperatures of molecular 
clouds, CO is frequently used as a tracer for H2 in extragalactic 
observations. The velocity-integrated CO line intensity (Wco) is 
converted to a molecular gas column density Ai'H2 using the conver- 
sion factor, Xco ■ Within the Milky Way and similar spiral galaxies, 
Xco appears largely consistent. In lower metallicity environments, 
such as dwarf and high-z galaxies, or in high-surface density envi- 
ronments, such as the Galactic Center and sub-millimeter galaxies, 
Xco is thought to deviate strongly from the standard value. 

Because of its importance in observations, Xco as a func- 
tion of environment is an active area of observational research 



(e.g.llsraellig??! : iLerov et al.l2009t iGratier et al.lbOIOt iLerov et al.l 



I2OII , [). Simulations of CO and H2 have also been used to 
asses the variation of Xco with the local gas properties. In a se- 
ries of papers inc l uding Glover et al. I teoioh .lGi over &Lowl(l20Ilh 
and IShettv et alj ( I2OI lal fo). Xco was studied using magneto- 
hydrodynamic simulations that followed the chemical network of 
CO and H2 within turbulent molecular clouds. Simulating CO 
requires high resolution and detailed radiative transfer but can 
be do ne assuming equilibrium JWolfire et al.ll201ol : iGlover & Lowl 
l20Ilh . Therefore, in galaxy scale simulations it i s often cal- 
culated during post-processing. In iFeldmann et al.l ( |20I2|) . the 
CO model from small-scale simulations was combined during 
post-p rocessing with simu l ations o f galaxies compu t ed wit h H2 
as in 'Gnedin & Kravtsov| ( |20I ih . iNaravanan et al.l ( 1201 ih and 
iNara vanan et al. (2012) combined detailed radiative transfer and 
equilibrium models of CO and H2 in the post-processing of iso- 
lated disk galaxies and mergers. By including a model for CO abun- 
dance in our simulations as done in these works, we would be better 
able to compare the molecular gas in our simulation to recent ob- 
servations of molecular gas in galaxies. Given the complexity of the 
the problem and the necessity of invoking an additional model for 
CO during post-processing, we have left the study of CO to future 
work. 



5 CONCLUSION 

In this work, we described a method for calculating the non- 
equilibrium abundance of H2 in SPH cosmological simulations. 
Using this method, we were able to reproduce the observed tran- 
sition from atomic to molecular hydrogen as a function of column 
density. We also showed that dependence of H2 abundance on vol- 
ume density and metallicity is consistent with previous studies. 
We demonstrated that this method produces consistent results for 
gas particle masses of approximately 10'' Mq and less. The inclu- 
sion of H2 and shielding enabled gas to reach lower temperatures, 
which resulted in us being able to simulate the cold phase of the 
ISM. In addition to changes to the themodynamics of the gas, the 
inclusion of H2 allowed us develop a more physically motivated 



star formation recipe in which star formation is directly linked to 
the local H2 abundance. By making star formation dependent on 
H2, star formation also became dependent on the local metallicity 
and LW radiation. 

Using this new method, we created the first cosmological sim- 
ulation with H2 of a galaxy integrated to a redshift of zero. This 
dwarf galaxy with H2 was compared to a galaxy produced with- 
out H2, in order to study the effect the inclusion of H2 has on the 
stellar and gaseous content. We found that when H2 was included 
in the simulation, a larger baryonic mass was in the disk of the 
galaxy at a redshift of zero, resulting in a brighter and more gas 
rich galaxy. The inclusion of H2 also resulted in more star forma- 
tion at later times and a bluer galaxy. Finally, the simulation with 
H2 had a larger distribution of star formation. 

These changes can all be related back to changes to the ISM 
temperature and density distributions. The presence of shielded gas 
and, to a lesser extent, the additional H2 cooling resulted in the 
increased formation of cold gas. It also resulted in a clumpier ISM 
and allowed regions dense enough for star formation to form at 
larger radii. 

These simulations demonstrate the importance of modeling 
H2. Changes to the thermodynamics of the ISM have a strong ef- 
fect on the mass content and on the spatial and temporal distribution 
of star formation. The spatial distribution of stars is additionally im- 
portant because it can affect both the baryonic and DM distribution 
of mass through feedback. 

In future work, we will examine how the interaction between 
feedback, star formation, and H2 affect disk properties and star 
formation in different masses of galaxies and at different redshifts. 
While they are sensitive laboratories for star formation, dwarf 
galaxies have only small amounts of H2. Larger, more metal rich 
galaxies contain more H2 and so simulations of them could be 
more affected by the larger fraction of cold, clumpy gas. We, there- 
fore, intend to study the effect H2 on star formation and feedback 
over a range of galactic masses. Furthermore, a more in-depth ex- 
amination simulations of galaxies with H2 at other redshifts and 
across their evolution will give insight as to how the evolution of 
the gas metallicities and surface densities and the accretion of mat- 
ter affect the presence of H2 and H2-based star formation. 
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APPENDIX A: ISOLATED SIMULATIONS 

We used an isolated Milky Way-like disk galaxy (virial mass equal 
to 10^^ M0) to test our implementation of H2 formation and de- 
struction at varying resolutions and metallicities and to test our esti- 
mate of the local LW flux. Isolated simulations have the advantage 
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that they can be simulated at much higher mass resolution than cos- 
mological simulations of similar mass galaxies. We choose an ap- 
proximately Milky Way-mass simulation in order to compare our 
results to observations of the Milky Way. 

In order to test our H2 implementation, we created the initial 
conditions for a set of Milky Way-like disk galaxies of different 
mass-resolutions and metallicities in the following manner. Prior to 
the addition of H2, we allowed a stabl e disk galaxy to form fr om a 
DM halo con taining a gas cloud, as in lKaufmarm et al] t2007l) and 
IStinson et all llooe). To do this, we beg an from a live, equilib rium 
Navarro-Frenk- White (NFW) DM halo dNavarro et al.lll997l) con- 
taining a cloud of gas. The cloud of gas initially shared the same 
NFW profile as the DM, had a uniform rotational velocity, and had 
a temperature distribution such that it was in hydrostatic equilib- 
rium prior to cooling. We then integrated the simulation for 1.0 
Gyr without H2 while allowing the gas cloud to collapse through 
ionization and metal-line cooling and to form stars according to our 
non-H2 star formation recipe. After 1.0 Gyr, a stable galaxy with 
a stellar and gaseous disk had formed inside the DM halo. We ad- 
justed the metallicity of the gas to reflect the range of metallicities 
we were interested in examining by adding a constant metallicity 
offset to the disk gas. We chose the metallicity offset such that the 
mean metallicity of the gas in the disk of the galaxy equaled in turn 
1.0, 0.3, and 0.1 Z©. 

The set of disk galaxies of varying metallicities and mass res- 
olutions created using the above methodology became the initial 
condition for our simulations used to test the H2 implementation. 
Table lAll contains a list of these galaxies and their parameters. For 
each of these galaxies, we enabled H2 in the galaxy and continued 
the simulations for a further 100 Myrs. This time period was chosen 
to be sufficiently larger than the star and H2 -formation time scales 
that the gas and star formation in the disk could reach equilibrium. 
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Table Al. Isolated Milky-Way like simulations used to test the code and their parameters. 
When calculating the mean smoothing length, {h), of the particles in the disk, we selected 
particles with densities of 0. 1 amu/cc or greater. The gravitational softening length was 
206 pc in all of the simulations. 



Particle Mass [Mq] 



Name Dark Star Gas (h) Mean Metallicity 

in Disk [pc] in Disk [Z/Zq ] 



MWk 10^ 4.0 X lO"' 1.4 X 10*' 278 1 

MWmr 10'^ 4.0 x 10"' 1.4 x 10^ 176 1 

MWmr0.3z " " " " 0.3 

MWmrO.lz " " " " 0.1 

MWhr 10^ 4.0 x 10^ 1.4 x 10* 82 1 
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